!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! the subroutine is used to calcaulate V, theta along DE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! coded by : B. G.
! created  : 2015-07-24
! revised  : 2015-08-05 revise error control method
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! include  :
!   subroutine : fta
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!   C2,C3   Lagrange mutiplier coefficients
!   thetad  flow(streamline) angle at point D
!   Md      Mach number at point D
!   yi~Vi   properties at prescribed point
!   theta0,Mstep  theta and Mach step
!   *temp*  temp variables
!   flagMg,flagMg0  flag if the current error greater than or equal to previous
!           one in M loop, 0-no, surrent theta loop and previous theta loop
!   mul*0,mul* the error product of previous and current loop, m-Mach, t-theta
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning:
!   1. angle in rad, not degree
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine Vtheta(C2,C3,thetad,Md,yi,Vi,thetai)
   use VariableDef
   implicit none
!   given vatiables
   real*8::C2,C3,thetad,Md,yi
!   temporary variables
   integer::ktemp,ktemp0,flagMg,flagMg0
   real*8::theta0,Mstep,thetatemp,Mtemp,alphatemp,temp1,temp2,roui,Mi
   real*8::mulm0,mulm,mult0,mult
!   calcaulate variables
   real*8::thetai,Vi
   
!   if (thetad==0) then
!      thetatemp=0.0
!      do Mtemp=Md,20,0.0005
!         Vi=Mtemp*sqrt(g*Rg*Tt/(1+(g-1)*Mtemp**2/2))
!         roui=Pt*(1+(g-1)*Mtemp**2/2)**(1/(1-g))/(Rg*Tt)
!         alphatemp=asin(1/Mtemp)
!         temp1=Vi*cos(thetatemp-alphatemp)/cos(alphatemp)
!         temp2=yi*roui*Vi**2*(sin(thetatemp))**2*tan(alphatemp)
!         if ((abs(temp1-C2) <= abs(0.0002*C2)) .and. (abs(temp2-C3) <= abs(0.001*C3))) return
!         if ((temp1 > C2) .or. (temp2 < C3)) exit
!      end do
!   end if
   
   theta0=-thetad/1e4
   Mstep=1e-4
   mult0=1000000
   ktemp0=0
   flagMg0=0

   do thetatemp=thetad+theta0,0,theta0   ! theta loop begin at thetad+theta0
      mulm0=1000000
      ktemp=0
      flagMg=0
      do Mtemp=Md+Mstep,20,Mstep   ! M loop begin at Md+Mstep
         Vi=Mtemp*sqrt(g*Rg*Tt/(1+(g-1)*Mtemp**2/2))
         roui=Pt*(1+(g-1)*Mtemp**2/2)**(1/(1-g))/(Rg*Tt)
         alphatemp=asin(1/Mtemp)
         temp1=Vi*cos(thetatemp-alphatemp)/cos(alphatemp)
         temp2=yi*roui*Vi**2*(sin(thetatemp))**2*tan(alphatemp)
         mulm=abs(temp1-C2)*abs(temp2-C3)
         ktemp=ktemp+1
         ! store minimum error
         if (mulm < mulm0) then
	         mulm0=mulm
	         mult=mulm
            flagMg=0
	      else
            flagMg=flagMg+1
         end if
         ! if only one crossed, next theta loop(temp1 > C2) .or. (temp2 < C3)
         if (temp1 > C2) exit   
      end do   ! M loop end
      ! compare with given value and return, begin
      if ((temp1 > C2) .and. (temp2 < C3)) then   ! if all of them crossed
         if (mult < mult0) then   ! if current error < previous one (theta loop), return this
            thetai=thetatemp   ! return current theta
            if (flagMg > 0) then   ! if current error >= previous one (M loop), return previous one
               Mi=Mtemp-flagMg*Mstep
               Vi=Mi*sqrt(g*Rg*Tt/(1+(g-1)*Mi**2/2))
               return
            else   ! if current error >= previous one (M loop), return this
               ktemp=ktemp+1
               ktemp=ktemp-1
               return
            end if
         else   ! if current error >= previous one (theta loop), return previous one
            thetai=thetatemp-theta0   ! return previous theta
            if (flagMg0 > 0) then   ! if current error >= previous one (M loop), return previous one
               Mi=Md+Mstep+(ktemp0-flagMg0)*Mstep
               Vi=Mi*sqrt(g*Rg*Tt/(1+(g-1)*Mi**2/2))
               return
            else   ! if current error >= previous one (M loop), return this
               Mi=Md+Mstep+ktemp0*Mstep
               Vi=Mi*sqrt(g*Rg*Tt/(1+(g-1)*Mi**2/2))
               return
            end if
         end if
      else
         mult0=mult
      end if
      ! compare with given value and return, begin
      ktemp0=ktemp
      flagMg0=flagMg
   end do   ! theta loop, end

end subroutine Vtheta
